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We present an analytical solution for the vortex lattice in a rapidly rotating trapped Bose-Einstein 
condensate (BEC) in the lowest Landau level and discuss deviations from the Thomas-Fermi density 
profile. This solution is exact in the limit of a large number of vortices and is obtained for the cases 
of circularly symmetric and narrow channel geometries. The latter is realized when the trapping 
frequencies in the plane perpendicular to the rotation axis are different from each other and the 
rotation frequency is equal to the smallest of them. This leads to the cancelation of the trapping 
potential in the direction of the weaker confinement and makes the system infinitely elongated in 
this direction. For this case we calculate the phase diagram as a function of the interaction strength 
and rotation frequency and identify the order of quantum phase transitions between the states with 
a different number of vortex rows. 

I. INTRODUCTION 

Rapidly rotating Bose-condensed gases constitute a novel class of many-body systems where the ground state 
properties are governed by a collective behavior of nucleated vorticesi*^. A harmonically trapped dilute Bose-Einstein 
condensate (BEC) strongly confined in the z direction, is essentially two-dimensional in the (x, y) plane. When the 
rotation frequency along the z axis becomes close to the trapping frequencies in the x and y directions, the BEC 
gas can be described as a system of interacting bosons in the lowest Landau level. The single-particle Hamiltonian 
is similar to that of a charged particle in a strong magnetic field, and the regime of fast rotation of neutral bosons 
presents an analogy with Quantum Hall Effect. Due to the presence of remaining harmonic trapping, the lowest 
Landau level (LLL) is not degenerate. However, analytic properties of the LLL wave functions generate an effective 
long-range interaction between the bosons, which results in an interesting physics. 

If the rotation frequency is not very close to the trap frequency, then the number of vortices is much smaller than 
the number of particles. Under these conditions the system is in the so-called mean-field Quantum Hall regime and 
can be described by a macroscopic wavefunction ^'(r) in the lowest Landau level. In this limit the vortices generically 
arrange themselves in a lattice. An increase in the rotation frequency increases the number of vortices and eventually 
it becomes comparable with the number of particles. This leads to melting of the vortex lattice and to the appearance 
of strongly correlated state^i-i^. The "mean-field Quantum Hall regime" for trapped bosons has been introduced by 
Ho^ and studied in a number of papers where the vortex lattice structures have been obtained numerically in the case 
of a circularly symmetric trapping potcntial'*'^-=i£i2.. 

In this paper we consider a rotating BEC in the lowest Landau level in the mean-field regime and obtain an analytical 
solution for the vortex lattice of the harmonically trapped symmetric 2D gas. This solution is exact in the limit of 
a large number of vortices, and we discuss deviations from the Thomas-Fermi density profile. We then turn to the 
case of the "narrow channel" geometry, which is realized when the confining frequencies in the x and y directions are 
different, and the rotation frequency is equal to the smallest of them. Then, in the rotating frame, the gas becomes 
extremely elongated in the direction of the smaller frequency, as has been demonstrated in the ENS experiment 
with thermal bosons^. This is an extreme case of a rapidly rotating 2D gas in an asymmetric harmonic potential, 
discussed in relation to the density profile of the gas and the density of vortices in Refi^. Some vortex structures of 
the asymmetric rapidly rotating BEC have been discussed and calculated in Refsj ^^'^^i^^'^"^ . In the narrow channel 
geometry, the excitation spectrum of a weakly interacting BEC without vortices exhibits a "roton-maxon" structurei^. 
The phase transition to the state with a vortex row occurs when the roton energy reaches zero under an increase in 
the rotation frequency or in the strength of interaction between the bosons. A further increase of these quantities 
increases the number of vortex rows through a set of quantum phase transitionsi^ii^. We classify these transitions and 
find an analytical solution for the vortex lattice in the narrow channel, which is exact in the limit of a large number 
of vortex rows. 
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II. GROSS-PITAEVSKII EQUATION IN THE LOWEST LANDAU LEVEL. SOLUTION FOR A 

SYMMETRIC HARMONIC POTENTIAL 



Consider a system of bosonic neutral atoms strongly confined in the z direction by an external trapping potential 
with frequency such that the bosons are in the ground state of the harmonic well and become essentially 
two-dimensional in the (x, ?/)-plane. The bosons are confined in this plane by a harmonic trapping potential ^(r), 
with r = {a;,y}, and the trap is rotating around the z axis with frequency 17. In the mean-field Quantum Hall limit, 
we assume to zero order that all particles are in the same macroscopic quantum state described by the wavefunction 
■!/'(r). In the rotating frame the Gross-Pitacvskii equation for '^(v) reads: 

|^V + <?l^lV + V^(r)^-f74^ = MV', (1) 
2 m 

where p is the momentum operator, m is the particle mass, is the operator of the orbital angular momentum, 
/i is the chemical potential, and is normalized to the total number of particles iV. Equation ([1]) is obtained for a 
short-range interaction between particles, and the 2D coupling constant g can be expressed through the 3D scattering 
length flg. If the harmonic oscillator length in the z direction, 1^ = ^Jh/miUz, is much larger than \as\ and the 
characteristic radius of intcrparticle interaction, then we have^^: 

ff = • (2) 

We will study Eq. ([T|) projected onto the lowest Landau level. A general procedure of obtaining the projected 
equation is described in the Appendix, and here we outline the method. 

The single-particle Hamiltonian for rotating neutral atoms is equivalent to the Hamiltonian of a charged particle 
in a uniform magnetic field B along the z axis. The field is such that half the cyclotron frequency ojc = -B/2m (in 
units of charge divided by the light velocity) is identified with the rotation frequency fi, and the vector-potential 
in the symmetric gauge is A = [B x r]/2 = r7i[r2 x r]. In the case of a symmetric external harmonic potential 
V{r) = niLiP'r^ 12 the single-particle Hamiltonian reads: 

li = ^-^i. + \mu?r'' = ^(P - eA/c)2 + )^m[u? - V?y , (3) 

At the critical rotation frequency 17 = w the residual confining potential vanishes, and the harmonic oscillator length 
I — {h/mujY/'^ of the initial trapping potential V{r) coincides with the "magnetic length" (/t/mfi)^/^. One then has 
an " infinite plane" geometry actively studied with respect to the ground state of interacting boson o^'^'i^^'^^i^"'^^ . 

Below the critical rotation frequency, Q. < oj and — fi) <C 17, the energy eigenstates are associated with the 
Landau levels of a charged particle in the uniform magnetic field, and the presence of the residual confining potential 
lifts the LLL degeneracy. A complete set of LLL eigenfunctions is given by 

4'„(z,z) = -^=exp [-^1 , (4) 
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with z,z = (a; ± iy)/l, and n being a non-negative integer. An arbitrary function in the LLL can be written as a 
linear superposition of the LLL eigenstates and represented in the form: 



/(^) 



I 



zz 
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^'(z,z) = — -^exp , (5) 



where f{z) is an analytic function of z. The projection operator onto the LLL is written as 

P = J2\n){nl (6) 

n 

where {z\n) — \E'„(z,z) is given by Eq. ([4]). Acting with the operator P on an arbitrary function (j){z,z) = 
[f{z, z)/l] exp(— zz/2) one obtains: 

z)^Y. *"(^' ^')<t>{z\ z')dz'dz' = [f{z)/l] exp(-zz/2), (7) 
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with dz'dz' = dx' dy' /I , and an analytic function /(z) which is given by 

f{z) = ^ J dz'dz' f{z', z') exp(zz' - z'z'). (8) 

This formalism was introduced by BargmannSS. and developed by Girvin and Jach^ in relation to Quantum Hall 
physics. 

In the case of interacting bosons one can still project the many-body Haniiltonian onto the lowest Landau level 
provided that the interactions are much smaller than the cyclotron gap, i.e. n2Dg ^ 2hfl, where n2D is the two- 
dimensional particle densityi^. Acting with the LLL projector P onto the Gross-Pitaevskii equation ([T]) results in 
the projected equation (see RefSf^iiS and Appendix): 

h{uj - n)zdj{z) + dz'dz'\f{z')\^f{z') exp(zz' ~ 2z'z') = fif{z), (9) 

where 11 = 11 — Huj, and the function [f{z)/l] exjp{—zz/2) is normalized to unity. Equation ([9]) has a simple solution 

Uz) = (10) 
V 7rn! 

which corresponds to the chemical potential and energy per particle given by 

Ng (2n)\ 

^^ = n{u:-n)n+—j-^, (11) 
^" o\ _u ^9 (2n)! 

- = n{u-n)n+—^j-^. (12) 

For n = 0, equation describes the ground state without vortices, and for n ^ it gives excited states with a 
(multicharged) vortex at the origin. The LLL approximation is valid when En ^ hiuN . The spectrum En (fT2)) has 
a roton shape with a local minimum at a certain value of n. In the limit of large n we have n\ w \/27m(n/e)", and 
Eq. (fT2)) is reduced to 

— = h{u:-^l)n+—^. (13) 

The local energy minimum is obtained for 

n = rtn. = 

47r V ^(w - ri) 



n = nQ = — I — — , (14) 



and from Eq. (|13p we find 



(15) 

The giant vortex state at n = no is a metastable state and it can have a relatively long lifetime. One can think of 
creating this state in dynamical studies and identifying it through the measurement of the density profile^. 

For Ng/l'^ ^ h{uj — Vl) the ground state represents a vortex lattice. At the critical rotation frequency = one 
finds an exact solution describing this lattice. Consider the function 

/o(z) = ^^^^1 (W&i, r) exp(^zV2«,), (16) 
where S is the surface area, r — u ^ iv^ and i9i is the Jacobi theta-function^-^ given by 

^ oo 

MC,t) = - y (-l)"exp{z7rr(n + 1/2)2 + 2<(n + 1/2)}. (^^^ 

n— — oo 

The Jacobi theta-functions are analytic in the complex plane and have zeros at the points = mr + tottt, where n, m 
are integers. The function fo{z) vanishes at the lattice sites nbi + mb2, with 62 = 6it. These points correspond to 
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the vortex locations. The parameter bi can always be chosen real so that the area of the unit cell is Vc = b^v. The 
absolute value of /o(z) is polar symmetric 



|/o(z)| exp 



2vcP 



(18) 



For any lattice with a fixed elementary cell area one has Vc — tt, and the normalization coefficient in Eq. ()16|) is chosen 
such that the function [fo{z)/l] exp{—zz/2) is normalized to unity. The function fo{z) is an exact solution of Eq. ^ 
for Lu = fl. It has a constant envelope and a periodic vortex structure. The minimum energy is obtained for the 
triangular lattice, where r = exp27ri/3, v — V3/2, and bf = 27r/-\/3. The chemical potential is then given by 

aNg 



with a ^ {3^^^/2) Y,b ^(-1)"'' exp{-7r2(62 + c2)/46?} = 1.1596, and b = 2m,c = 2p being even integers. 
In the case of > il, a general solution of equation © can be represented as 



f{z)^{2vf'^ E (-irff(«)9'^' exp 



ITT Z 



(19) 



where a = 2n + 1 are odd integers, q = exp[i7rr/4], and g{a) is a differential operator acting on a. For g{a) = 1 
one recovers the solution (fTB)) . Substituting the trial function (fT5)) into equation ^ for a triangular-like lattice (the 
lattice that becomes exactly triangular for a; = f2) we obtain: 



{fi* - A+A)g{a) -=- E 3(0 - b)gia - c)g{a - b - c) exp 



-451 + 



(-!)• 



xg" exp 



-az 



0. 



(20) 



Here (3 — Ngl{Ph{u) — fi)), /i* = ji/h{u) — 57), and we introduced the operators 

A,A+ = —a±^ — . 

2bi TT oa 

For large f3 we have /i* ^ 1 and an approximate solution for g{a), which describes the vortex structure with a high 
accuracy, turns out to be 



5(a) 



1 



i?2 _ A+A e(i?2 - it^)^ 



(21) 



where Q is the Heaviside theta function, R = V/i*, and we will see below that it is the radius of the condensate cloud 
in imits of /. Substituting the solution ([2T|) into Eq. (fT9| we obtain after some algebra: 



{2v)^ 



n= — 00 /c=0 



2fe/2/,r 



(2n + l) exp{-7rw(2n+ 1)2/4}, (22) 



where Hk(w) are Hermite polynomials. 

Equation ((2T|) is obtained taking into account that the leading contribution to the sum over b and c in Eq. ([20]) 
comes from small values of m and p, since already the contributions of terms with |m| > 2 or |p| > 2 are exponentially 
small. Provided that the dependence g{a) is smooth, which is the case for large i?, we may consider large a and omit 
b and c in the arguments of the g operators in Eq. ([^0]) . This immediately gives Eq. (PT|) . From the condition that 
the function [f[z)/l] exp(— zz) is normalized to unity we find a relation R — {2aj3/iTY^^, in agreement with Rcfs.^'^^, 

For the angular-averaged particle density, i.e. the density averaged over the azimuthal angle (p [z = r exp iip) we 
then have: 



n{r) j |^(z, z-)pg = nM^vf" E E " ^) ("1) 



{[n(n-l)+m(m-l)]/2]} 



J2k 



2'=(fc 



M2 



Hk J — {2n + l) Hk W— (2m + l) exp{-^t;(2n + 1)74 - 7r?;(2TO + 1)74} exp(-r2), 



(23) 
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FIG. 1: Angular-averaged density n{r) in units of n2D versus r (in units of I) for R — 11. The solid curve shows the result of 
Eq. 1231 and filled circles the results of the variational calculation (see text) . 
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FIG. 2: Condensate wave-function \ip{x,y)\'^ for R = 11. Coordinates x (horizontal line) and y (vertical line) are given in units 
of 

where n2D = 2N/ttPE? is a characteristic 2D density in the central part of the cloud. The angular-averaged density 
calculated from Eq. (|23p for i? = 11 is shown in Fig. [TJ In the entire region of r it coincides with the numerical 
result obtained by expanding the condensate wavefunction in terms of the single-particle LLL states ([4]) and using 
a variational approach for finding the coefRcients of the expansion. This demonstrates a very high accuracy of our 
analytic solution. The structure of the vortex lattice for i? = 11 is shown in Fig. [51 

The angular- averaged density represents oscillations on a length scale of the magnetic length Z, on top of a slowly 
varying envelope. Averaging the density over the oscillations, that is averaging n(r) (or just the density n{r)) over a 
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FIG. 3: The coarse grained density in units of n2o versus r, calculated from Eq. (|26p for R — 11. The dashed curve shows the 
Thomas- Fermi inverted-parabola shape (|27|l . and r is given in units of I. 



distance scale much larger than I, gives the coarse grained density: 



ncg{x,y) — — j dSx ( d5yn{x + 5x,y + 5y) (x\f{z + 5z)\^ ex.-p[-~{z + 5z)[z + 5z)], 
^ Jo Jq 



where R C ^ 1. Using the function f{z) of Eq. (fTO|) we have: 



(24) 



ncg{x,y)(x— / dSx dSy ^ (-l)"+™ff(a„)3(a„) 
'-Jo Jo „ 



n.m— — oo 



exp 



iTTToi iTrT*a'L in, ^ , tt , ^ x „/ c \9 

—r^ -7-^ + -ri^ + 5x){an -am)--r{y + 5y){an + a™) - 2{y + 5yf 

4 4 Oi Oi 



(25) 



with a„ — {2n + 1), and a„i — (2m + 1). Integration over dSx gives n = m and transforms the y-dependent part of 
Eq. (USD to 



1 r °° 



exp ■ 



n— — oo 



/TT/ r ^ 7r(2n + 1) 
V2(, + .,) + ^^ 



One then clearly sees that the integration over dSy/C is equivalent to replacing the summation over n by integration. 
Thus, in order to obtain the coarse grained density from Eq. we have to put n = m and integrate over dn. This 
yields: 



fc=0 



i?2 ; r(i?2) 



"R2r(i?2) 



(26) 



where T{R^) and r(i?2^r2) are the Gamma function and Incomplete Gamma function, respectively. 

For r < R and (i? — r) 3> 1 we may put r{R^, r^) ~ r(i?2) in the right hand side of Eq. and omit the second 
term which gives a correction of the order of \/ R or smaller. This leads to the expected Tomas-Fermi density profile: 



(27) 



ncg{r)^n2D [l--.] Q{R-r) 



i?2 



For r > R and (r — i?) 1 we have at large R: 



r(i?2,r2) = 



1 



r2 - i?2 (r2 _ i?2)J 



(r2)^" exp(-r2). 



7 



Then, using an asymptotic expression r(i?^) = ^JQjtJB? {R^)^^ exp(— i?^), we obtain that the density decays expo- 
nentially: 

and is practically zero for r > The coarse grained density versus r at i? = 11 is displayed in Fig. O At the 

Thomas- Fermi border, r — R,we have ficg = n2D / \/2ttFP. The validity of the Thomas- Fermi inverted-parabola shape, 
in general, requires the inequality (i? — r) ^ 1. Nevertheless, in Fig. [3] we see that for i? = 11 the Thomas-Fermi 
formula works well already for r < 10. 

Deviations from the Thomas-Fermi density profile of ncg(r) have been studied in Ref.^ by using the variational 
procedure. Here we present an analytic solution and show that it describes very well the density profile, including the 
non-Thomas-Fermi part. 

III. NARROW CHANNEL GEOMETRY 

Let us now consider an anisotropic confining potential ^(r) = miuji'^x'^ + w^2/^)/2, with ojy < oj^- At the critical 
rotation frequency U, — ojy^ the centrifugal force cancels the confining potential in the y-direction. One then has a 
quasi-one-dimensional system in the rotating frame, which is usually refered to as the narrow channel geometry. The 
system is infinitely elongated in the y direction and is confined by a residual transverse potential 'rn(uj'^ — f2^)a;^/2 in 
the X directioni^. 

After the transformation to the Landau gauge, — s- 4* exp(imf2a;y//i), a complete set of eigenfunctions of non- 
interacting particles in the lowest Landau level of the narrow channel is given by 

where L is the length of the system in the y direction in units of I = {h/m^lY/"^ , = (w^ -I- 3ri^)/4, fc„ = 27rn/i 
with n being an integer, and we introduced dimensionlcss coordinates 

~ X _ 

2; = --^; Q = x + iy. (30) 

Thus, the wavefunction of any state in the LLL can be written in the form: 

vf = [/(O/f] exp . (31) 

The projection operator onto the LLL is given by Eq. (I6|), and acting with this operator on an arbitrary function 
C) — [/(C; C)/'] exp(— y^) wc obtain an analog of equations ([7]) and ([5]): 

P0(C,C) = E S / MUWn{C'X'm'X')dx'dv' = [}{0/l]eM-f): (32) 
/(C) = - / dC'dC7(C', CO exp{CC' + C'V2 - CC' - CV2}, (33) 



TT _ 

where dC'dC,' = {n/n)dx'dy' /P. 

The Gross-Pitavevskii equation projected onto the lowest Landau level in the narrow channel has the form (see 
Appendix for details): 

- fiu;o/"(C) + dC'dClfiCTfiC) exp (^~2CC + CC' + + ^ " y) - A/(C), (34) 

with ujQ ~ r2(f2^ — il^)/(2il^) ^ fl being proportional to the frequency of the remaining confinement in the x direction. 



/i = /i — ftw, and the condensate wavefunction [,f{C)/l] exp(— y^) being normalized to unity. In analogy with Eq. (|19p 
let us again search for the solution of the form 

(2z;)V4 ^ /. i2n + ll znCi2n+l)\ 
/(C) = ^ 2^ (-1) 5(2n + 1) exp l^iTTT + j , (35) 
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where bf = 2it/V3, v = a/S/Z and r = exp(27ri/3). Equation p5|) describes the structure with odd number of vortex 
rows, with the central row at y = 0. Using instead of in Eq. (|35p . which corresponds to the replacement 
{2n + 1) ^ 2n, we obtain structures with an even number of vortex rows. Substituting Eq. (|35p into Eq. (|34p yields 



3(a) 



3^/^V2^(3 g{a ~ b)g{a -c)g{a-b- c) exp 



(-1) 



(36) 



where 13 = Ng/2y/TrhiiJoLP, fi* = ji/fiWQ, a = 2n + 1, 6 2m and c — 2p are odd and even integers, respectively. 
As well as in the symmetric case, at large (3 we find an approximate solution for g{a), which describes the vortex 
structure with a high accuracy: 



ff(a) 



/ 1 



\2a^j3 



1/2 



— dAm/iTV- 
bi 



a2 e 



2R 



(37) 



where we put ji* = 4i?^, and it will be seen below that R is the half-size of the cloud in the x direction (in units of 
I). From the condition that the condensate wavefunction is normalized to unity we obtain: 



(3a\/^/3/8)i/^ 



(38) 



Equation (|37p is obtained by putting the arguments of the g functions equal to a in the sum over 6, c in Eq. (j36p . 
Similarly to the symmetric case, the contribution of terms with high b and c in this sum (|m| > 2 or \p\ > 2) is very 
small, except for n very close to the border value biR/i:. The relative contribution of such n to the sum in Eq. (j35[) 
decreases with increasing R. Thus, the solution f{Q ([55]) with g{2n + 1) of equation ([57)1 becomes exact in the limit 
of large f3. The structure of the vortex lattice for (3 — 50 is displayed in Fig. |4l For very large (3 the number of rows 
for a triangular-like lattice is approximately equal to 2i?/(V36i/2) cx f3^^^ (see next Section). Then, our results lead 
to the Thomas-Fermi density profile in the x direction for the major part of the cloud, as explained below. Using 
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FIG. 4: Density profile \'ip{x,y)\'^ for (3 ~ 50. Coordinates x and y are given in units of /, and IV"!^ in arbitrary units. 

equations ((35)) and (|37p we define the line-averaged density n{y), i.e. the density averaged over the direction y of 
vortex lines: 



n{y) = - I mx,y)\^dx^ ^ (—-(2^+1) 



L 



2a(3bl 



xexp<j-2[2/-f ^(2n+l)]2 



(39) 



with {2nmax + 1) ^ 2biR/'K. In Fig. [5] we compare the result of Eq. (I39p with the result obtained by expanding the 
condensate wavefunction in terms of the single-particle LLL states (j29p and using a variational approach for finding 
the coefficients of the expansion. The comparison shows a very high accuracy of the found analytic solution. 
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y 

FIG. 5: Line-averaged density n(j/) in units of n2D versus y = x/l for a condensate in the narrow channel for {3 — 900, 
R = 8.8397. The sohd curve shows the result of Eq. (|39|l . the filled circles indicate the results of the variational calculation (see 
text), and the dashed curve the Thomas- Fermi inverted-parabola density profile. 



The line-averaged density n(jj) shows oscillations on a length scale I, on top of a slowly varying envelope. 
Averaging the density over the oscillations, i.e. averaging n{y) over a distance scale much larger than I, gives the 
coarse grained density. The averaging procedure is equivalent to replacing the summation over n in Eq. ([39]) by 
integration, and we obtain the following expression for the coarse grained density: 



ncg{y) = «2D 



'2 1 



2 1 



R 



dw(R^ - w^) exp{-2(y + wf} 

{erf[V2(^ + R)] - erf[V2(y - R)]] 



■exp{-2(y-7?)2} 



(40) 

where n2D — Suid/ARI is a characteristic 2D density, and = N/Ll is the one-dimensional density in the narrow 
channel. The coarse grained density versus y for R = 8.9 is displayed in Fig.[6l For {R— \y\) ^ 1, Eq. (|40|) immediately 
gives the expected Thomas-Fermi density profile: 



ncgiy) = n2D 1 



i?2 



(41) 



and for y — R we have ricg = 1/V ^ttR"^. As well as in the symmetric case, the inverted-parabola formula already 
works well not very far from the Thomas-Fermi boarder. In Fig.[n]one sees that for R — 8.9 this is the case for \y\ < 8. 
If \y\ > R and {\y\ — i?) ^ 1, then the density decays exponentially: 



ricgiy) = n2D 



27ri?2 4:{\y\~R)- 



.exp{-2(|y|-i?)2}. 



(42) 



The melting of the vortex lattice occurs when the number of vortices Ny becomes of the order of the number of atoms. 
The number of vortex rows at large /3 increases as and the spacing between the vortices is ~ I. Thus, the number 
of vortices is ~ P^^^L, and the melting transition occurs at the one-dimensional density nuj = {N/ LI) ~ (3^/^/1. 
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FIG. 6: The coarse grained density in units of n2D versus y 
shows the Thomas-Fermi shape. 



x/l, calculated from Eq. (|40p for R = 8.9. The dashed curve 



IV. PHASE DIAGRAM FOR A CONDENSATE IN THE NARROW CHANNEL 



In this Section we calculate the phase diagram for a rapidly rotating Bose-Einstein condensate in the narrow channel 
geometry. The phase diagram is obtained by numerical minimization of the energy functional 



£'//kJo = ^ fc^la/cP + /3 ^ afc+qa^,_qafc'afcexp 

k k,k',q 



(43) 



where we impose periodic boundary conditions along the y-axis and omit the index ri for momenta A;„. The energy 
is measured in units of Hluq and depends only on a single dimensionless parameter f3. The functional (|43p is obtained 
by substituting the condensate wavefunction 

i!{x,y) =J2(^k'i^k(x,y), (44) 

with '^k{S:,y) given by Eq. (|29|1. into the Hamiltonian for interacting bosons in the Landau gauge and integrating 
over dx and dy—. 

The coefficients Ofc were calculated by minimizing _E (|43p using a simulated annealing algorithm^^-. In general, there 
is an infinite number of coefficients ak in the variational wavefunction (I44|) . but the ones corresponding to large 
momenta are strongly suppressed due to the presence of the "kinetic energy" term in the energy functional. The 
normalization condition for the condensate wavefunction leads to the constraint \af. = 1. 

At /3 = the energy is minimized by setting all coefficients Uk with k equal to zero and gq ~ 1. This corresponds 
to the condensate density profile shown on Fig. [7^, which is a Gaussian with the half- width / in the x direction and is 
uniform along the y axis. This state remains the ground state for (3 < 4.9, and for (3 = 4.9 it transforms via a second 
order quantum phase transition into the state displayed in Fig.[7)D. In this state two extra components, ko and — fco, 
develop and the ordering wavevector fco has the value fco — 2.25 for (3 ~ 5.2. The three main components of this state 
are accompanied by nonzero, but much smaller components with higher fc which are multiples of fco. The critical 
value of f3 for this phase transition can be obtained analytically by minimizing the energy of the three-component 
wavefunction (|44)) . with fc = and fc — ifco^^. In this case the emerging state is seen as two rows of vortices^^, 
although including small components with higher fc it becomes a sort of corrugated state and can also be identified 
as a density wave. 

At /3 « 5.4, there is a first order phase transition from the density-wave state b) into the state with one vortex row 
(Fig. [T):). In this state the central component vanishes (ao = 0) and the wavefunction is characterized by two main 
non-zero components a±ko ~ 1/a/2- At /3 ~ 10 the ordering wavevector is equal to fco — 1.51. The energy of the 
purely two-component state is larger by a small amount than the energy of the state c) , which is especially visible 
near the phase transitions. 
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a) 



b) 



c) 



d) 



FIG. 7: Condensate wave-function \ip{x,y)\^ for different values of the interaction strength: a) /3 = 0, b) /3 = 5.2, c) /3 = 10, d) 
/3 = 19.2. Coordinates x (horizontal line) and y (vertical line) are given in units of I. 




I 
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e) f) g) 

FIG. 8: The same as in Fig. [7] for: e) = 30, f) /3 = 50, g) /3 = 100. 



Close to /? = 18.9, the state c) transforms into the state shown in Fig.[7tl and representing a density wave of vortices. 
This looks like the first order transition (see Ref.^^). However, the state d) becomes the ground state only for /3 > 18.88, 
whereas the state c) is the ground state for /3 < 18.85. In the narrow interval 18.85 < /3 < 18.88 our calculations yield 
a dynamically unstable corrugated state d), which is signaled by a negative sign of the compressibility. This is likely 
to mean that in this narrow range of (3 the system undergoes the phase separation. 

The state d) has five main components, including the central component at A: = and has the ordering wavevector 
which is equal to ko ~ 1.74 at ~ 19.2. This state turns into the state with two vortex rows (Fig. [5^) via a first-order 
phase transition at (3 ~ 19.8. 

For larger (3, there are phase transitions at /3 ~ 47.7 to the state with three vortex rows (Fig. [8f), and at /3 ~ 94 to 
the state with four vortex rows (Fig. [8|;). These transitions seem to be of the first order. The physical explanation 
for the possible absence of intermediate corrugated/density- wave states near the first-order phase transitions into the 
states with a larger number of vortices could be that starting from two vortex rows the system becomes rigid to 
corrugations in the transverse direction. So, increasing f3 we observe an increase in the number of vortex rows through 
first order transitions. For /? ~ 164 there is a transition to the state with five vortex rows, for (3 ~ 260 to the state 
with six vortex rows, and for /3 ~ 385 to the state with seven rows, etc. (see Fig IH] and Fig. [TT]) . Already for the 
state with 8 vortex rows, which emerges at /3 ~ 550, the composition of the rows looks like a triangular lattice (see 
Fig. [TUl) . For a large number j of the vortex rows, the Thomas- Fermi size of the cloud in the x direction, 27?, satisfies 
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h) k) 

FIG. 9: The same as in Fig.[7]for: h) /3 = 200, k) /3 = 300. 




-10 -5 5 10 



FIG. 10: The same as in Fig.[l]for: /3 = 600 (see text). 

the asymptotic relation ((38)) and is proportional to (3^^^. It is equal to the distance between the rows multiplied 
by (j + 1). Thus, the value of [3 corresponding to the transition from [j — 1) to j vortex rows obeys the relation 
Pj = f3j-i[{j + It works with a high accuracy, which is better than 0.5% for j > 8. 

The narrow channel geometry for rotating Bose gases was first considered in Ref;^^, where the roton-maxon structure 
of the excitations of the BEC without vortices has been found, and the phase diagram was presented with an emphasis 
on the first two transitions which can be calculated analytically. A numerical study of the related problem was done 
in RefplS., where corrugated states were discussed. However the analysis of quantum transitions in Ref.^^ stops at 
the appearance of the state with two vortex rows, although the states with up to 4 rows of vortices have also been 
observed. Here, we present a complete phase diagram and identify the nature of quantum phase transitions. The 
chemical potential as a function of (3 for /3 < 50 is shown in Fig. ll2[ indicating three first order transitions, one second 
order transition, and the above described peculiar transition between the states c) and d). 

The extension of the condensate wavefunction in the x-direction increases with increasing the interaction strength, 
and the two-dimensional density decreases. This decreases the average filling factor defined asv — N/Ny. The Gross- 
Pitaevskii equation gives a good description in the limit of large filling factors and we expect our picture to break 
down at very large f3. Eventually, when the number of particles becomes of the order of the number of vortices iV„, 
the vortex lattice melts through the phase transition to a strongly correlated state. The limiting case of extremely 
large /3 corresponds to the Laughlin state with v = 1/2, and it was discussed for the narrow channel with periodic 
boundary conditions in the x direction in Ref.-^'. 
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FIG. 11: Zero temperature phase diagram for a rapidly rotating condensate in the narrow channel. Solid vertical lines indicate 
the points of first order transitions, and the dashed line the point of the second order transition. The bold solid line shows the 
transition between the states c) and d) (see text). The numbers from to 6 stand for the number of vortex rows in a given 
range of ^, and the filled areas correspond to corrugated/density- wave states. The letters from a) to k) indicate the figure in 
which a given vortex state is shown. 

15 



12.5 
10 



7.5 

3, 



2.5 




10 



20 



30 



40 



50 



/3 



FIG. 12: Chemical potential in units of fiwo as a function of /?. The dotted line indicates the transition between the states c) 

and d) (sec text). The insets show the dependence /i(/3) in the vicinity of the quantum transitions eX ji = 4.9 (upper inset) and 
at /3 ~ 18.9 (lower inset). The dashed lines in the insets indicate the derivative 9/i/9/3 in arbitrary units. 



V. CONCLUSION 



In conclusion, wc found an analytical solution for the vortex lattice in a rapidly rotating BEC in the LLL. This 
solution is asymptotically exact in the limit of a very large number of vortices, and we discuss non-Thomas-Fermi 
effects in the density profile. The results are obtained for two limiting cases, circularly symmetric BEC and narrow 
channel geometry. In the latter case we present a complete phase diagram and identify the order of quantum phase 
transitions occurring under an increase in the interaction strength and/or rotation frequency and resulting in an 
increase in the number of vortex rows. 
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Appendix 

Let us give a detailed derivation of the projection of the Gross-Pitaevskii equation for a rapidly rotating Bose- 
condensed gas onto the lowest Landau level22.. The gas is confined in an asymmetric harmonic potential ^(r) — 
m{u)lx'^ + a;^y^)/2, and we assume without loss of generality that ujy < tOx- 

In the symmetric gauge the single-particle Hamiltonian is similar to that of equation ([3|): 

H^lip-[nX r])2 + i(c^2 _ ^2)^2 ^ 1 (^2 _ ^2^y2^ (45) 

where we put h = m = 1. Drawing an analogy with a charged particle in a uniform magnetic field B, the rotation 
frequency Q is identified with half the cyclotron frequency and putting the particle charge and light velocity equal 
to unity we have = ujc = B /2. In complex coordinates the Hamiltonian (j45p rewrites as 



H{z, z) = -2dd + ujc{zd - zd) + ^tj^^z + ^('^^ " + ^^)' (^6) 



where Lof = uj^ + Uy and can be rewritten as ^ uj^ + {u}^ + ujy)/2, with u)^ = cj^ — uj"^ and ujy — ujy — 



Introducing the frequencies uj^ ~ \/ uj'?. + f "^^"" )^ and the dimensionless parameter = . / ("f^-H'^c)("t +"<=) -(-j^g 
unnormalized ground state wavefunction is 

(z, z|*o) = e^^"^'^^^e-^(°^'+*^'\ (47) 



where a = ^p'^i^J^ - ^c), b = ^{uj^ + ujc)/p^, and i/|*o) = ^t\^o)- 

The lowest Landau level in the asymmetric well is obtained by redefining j'l') = l^'o^') and requiring the Hamiltonian 
H which acts on 5') to depend only on a single variable representing a linear combination of z and z: 



P 

The eigenvalue equation acting on |^) then reads: 



\z. (48) 



P 



{E - = + [ut - u}i)u^'. (49) 



Let us define a dimensionless variable u' = u)^ — lo^ u so that the eigenvalue equation (|49|) becomes 

[E - c^+)# = ^ {-^" + 2v!^') (50) 
This is a Hermite equation (cj^ > uj^') with eigenfimctions 5'„(m') — Hn(u') such that 

{Z,z\^n) ^ NnHn{u')^o{z,z) 

with eigenvalues 

En = n{uj^ — ujf^) + uj^ . 



Introducing the quantity /z: 



cosh ^ = < 

UJr 



sinh /z = 1 



/ 4- 




+ UJc) 





















so that 

and using the relation 

du'du' 



{'fo\z,z){z,z\^o) 



~u'u' taiih fi{u''^ -\-u' 



' ^sinh 2/i \/sinh 2/i 



) = cosh^ 



^ tanh_£i y 



4,;, 



the normahzation factor is found to be 



ujc cosh 

The projector onto the LLL of an asymmetric harmonic well, Plll ~ J2n>o \ '^n){^n\, is 

{zi,Zi\PlLl\z2,Z2) = '^''^K ^-l(a.;+b.i+^+.i.-i)e-l(a.-i+b.iW^2^2) 



TTWc cosh /X 



ra>0 



Using the relation 



one obtains 



2" 1"^ 2 



{zi,zi\Plll\z2,Z2) = ^-^e-5^-i"^"'i -^-^Jf^-^y 

TTUJc 



^-i(az^+bz^+ui + ziZi)^-i(azi+bzi+iU+Z2Z2) 
sinh^ fi{u'^ +U2)+sinh 2^u\u2 

Any state — Plll\'^) in the LLL is a linear combination of LLL eigenstates: 

oo 



A^LLl) = X] Cn(2, ^l^n) = f{u')e- 



n=0 

where /(w') is an analytic function. 

Consider now the Hamiltonian to which we add a scalar potential V{z^ z): 

H{z, z) ^ ~2dd + ujc{zd ~ zd) + ^^uj^zz + ^(w^ _ ujI){z^ + z^) + V{z, z) 

Projecting the eigenvalue equation H\^) = E\'i') onto the LLL amounts to {zi, zi\PlllHPll 
E{zi,zi\Plll\-^). This gives 

+ - r 

I ^^^^-^g-i(a2|+bz|+i/j+Z2Z2)g--|(az|+62|+cj+2222)g- sinh^ p,{u'j^+u'2)+sinh2fj.u[u'2 

Hiz2,Z2)fiu',) = Ef{u[). 
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Writing explicitly H(z2, z-i) and changing the integration variables to u'^, u'2, we have 

^i^^J^^/fg- sinh^ ^uf f gSinh^ /xU2^gSinh2/i(«'iii2-«2U2) 

TT J ^ ^ 

(j^^^i-nW^) + 2u'J'{u'2)) + cot f + ^("2, "2)/(4)) = Ef{u[). (55) 
Using the Bargman identity 



sinh 2fi 



one finally obtains 



{E~Lot)f{u') = ^ {~f"{u) + 2u'f{u')) + ^ ; e^'"^' '^""/(^O, 

where the notation : V{u' ^ ^i^j^^f^u') ^ means that in V{u' ,u') the variable w' has been replaced by the operator 
sinh 2^1 ^'^'^ normal ordering has been made. In the case of the Gross-Pitaevskii equation the scalar potential 
is replaced by the non-linear term 

g{^LLL\^LLL) ^ gj(u)i[u)e ^ T ; 

SO that Eq. ((55|l becomes 

{E - ^+)/K) = i^(-/"K) + 2</'K)) + 

^sin^g-sinhV^' Jd^d^i'2e2'^'"'^'''"2'+'^inhV^%2sinh2p (gg) 

Ea. (|56p is a general form of the non- linear Gross-Pitaevskii equation projected onto the LLL of an asymmetric 
harmonic trap. 

Let us now concentrate on the two cases of interest, circularly symmetric geometry and narrow channel geometry. 
In the symmetric geometry we have w^, = Wy = w = wt, and Lof — cu^ ^ cut — Wc, 

sinh/i — > — —{Cbx — ^y) — * 0, cosh/i — > 1 

— 2iu' — *■ \/u}t — bJcPZ 

I ujuTc 
2. /- — — 00 

y {uJt - Ulc){<^x - Wy) 

SO that 



LOOJc 

U ^ I. j — —z — !■ lOO. 



Changing variables, u' z, equation reduces to 

{ut - c^,)zi/'(zi) + g'^ j dz2dz2e^^*('-''/^-'-''^f{z2).fiz2) ^{E- Ut)f(zi). (57) 
Finally, when w — > (critical rotation), i.e. uJt Wc, (f57|) becomes 

g'± j dz2dz2e^'^^^^''''/^-''''^fiz2)fiz2) = {E- LUc)f{zi). (58) 

Putting UJt =uj, ujc = ^, (E - UJt) = fx, and / -> {VN/l)f in Eq. we obtain Eq. @. 
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In the narrow channel geometry ujy — s- 0, and — lo^ 



sinh p, = cosh p 



-2iu' 

\ ~ 



where uj[ — ^ uj'^ -\- Changing variables, u' ^ u — z — z/ , equation (j56p reduces to 

where lo'I = ajf C^*^^" )^. Turning to the variable C = i-\/w"w2 , putting = i^, and noticing that = f2 where Cl 
was defined in Section III, we have 2a;" (w^ — a;c)/('^f +^c) — with wq introduced in Eq. ([34]). Then, after rescaling 
the function / as / ^ (VN /l)f, equation ([55]) transforms into Eq. ([M]) . 



a;a;(Iij^/(2wc 
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